#include <boost/algorithm/string.hpp>
#include <iostream>

#include "blosum.h"

using namespace boost;

Blosum::Blosum () {

}

void Blosum::setAlphabet(string alphabet  = "AGCT") {
    alphabet = to_upper_copy ( alphabet );
    int index = 0;
    for ( string::const_iterator c = alphabet.begin(); c != alphabet.end(); ++c ) {
        alphabetSet_.insert ( *c );
        alphabetMap_.insert ( make_pair ( *c, index++ ) );
    }
    alphabetLen_ = alphabetSet_.size();
}

void Blosum::fillSequencesFromFile ( const char* file_path ) {
	ifstream dataFile;
	string textLine;
	int lineCount=0;
	vector<string> container;

	dataFile.open(file_path);

   if(dataFile.is_open()) {	
		while(!dataFile.eof())
    	{
  		 getline(dataFile, textLine);
		 
		 if(lineCount == 0)
		{
			setAlphabet(textLine);		  
			cout<<"Wczytywanie z pliku..."<<endl;
			cout<<"alfabet: "<<textLine<<endl;
		}
		 else
		 {
			if (textLine!="") {
			 container.push_back(textLine);
			 cout<<"sekwencja nr: "<<lineCount<<" - "<<textLine<<endl;			
			}
		 }
		 lineCount++;
		}
  dataFile.close();
  setSequences(container);
  }
  else
    cout<<"Błąd: Nie można wczytać danych z pliku"<<endl;

}

void Blosum::setSequences ( const vector< string >& sequences ) {
    int length = 0;
    sequences_.clear();
    sequences_.reserve ( sequences.size() );
    for ( vector<string>::const_iterator i = sequences.begin(); i != sequences.end(); ++i ) {
        string sequence = to_upper_copy ( *i );
        for ( string::const_iterator c = sequence.begin(); c != sequence.end(); ++c ) {
            if ( !isInAlphabet ( *c ) ) {
                break;
            }
        }
        if ( length == 0 ) {
            length = sequence.size();
        } else if ( length != sequence.size() ) {
            break;
        }
        sequences_.push_back ( sequence );
    }
    sequenceLenght_ = length;
}

void Blosum::calculate() {

    vector< vector<int> > counts ( sequenceLenght_, vector<int> ( alphabetLen_, 0 ) );
    calculateCounts ( counts );

    matrix<int> step1Matrix ( alphabetLen_, alphabetLen_ );

	for (unsigned i = 0; i < step1Matrix.size1 (); ++ i)
        for (unsigned j = 0; j < step1Matrix.size2 (); ++ j)
           step1Matrix (i, j) = 0;

    int sum = calculateStep1Matrix ( step1Matrix, counts );


    matrix<double> probabilitiesMatrix ( alphabetLen_, alphabetLen_ );

	//zerowanie	macierzy prawdopodobienstw
	for (unsigned i = 0; i < probabilitiesMatrix.size1 (); ++ i)
        for (unsigned j = 0; j < probabilitiesMatrix.size2 (); ++ j)
           probabilitiesMatrix(i, j) = 0;

    vector<double> marginalProbabilities ( alphabetLen_, 0 );
    calculateProbabilitiesMatrix ( probabilitiesMatrix, marginalProbabilities, step1Matrix, sum );
    calculateBlosum ( probabilitiesMatrix, marginalProbabilities );

    testShowCounts ( counts );
    testShowStep1 ( step1Matrix, sum );
    testShowProbabilities ( probabilitiesMatrix, marginalProbabilities );
    testShowBlosum();
}

matrix< int > Blosum::getBlosumMatrix() const {
    return blosum_;
}

void Blosum::saveBlosumMatrix ( string file_path ) const {
	ofstream dataFile;
	const char *Cfile_path = file_path.c_str(); 
	string textLine;	
	dataFile.open(Cfile_path);

  if(dataFile.is_open()) {
   dataFile << "Macierz BLOSUM: " << endl;
   dataFile << "\t";
    for ( int i = 0; i < alphabetLen_; i++ ) {
        dataFile << indexToChar ( i ) << "\t";
    }
    dataFile << endl;
    for ( int i = 0; i < alphabetLen_; i++ ) {
        dataFile << indexToChar ( i ) << "\t";
        for ( int j = 0; j < alphabetLen_; j++ ) {
            dataFile << blosum_ ( i, j ) << "\t";
        }
        dataFile << endl;
    }
 dataFile.close();
 } else
  cout<<"Błąd: Nie można zapisać pliku z macierzą blosum!"<<endl;
}

void Blosum::calculateCounts ( vector< vector<int> > & counts ) const {

    for ( vector<string>::const_iterator sequence = sequences_.begin(); sequence != sequences_.end(); ++sequence ) {
        const string & seq = *sequence;
        for ( int i = 0; i < sequenceLenght_; ++i ) {
            ++counts.at ( i ).at ( charToIndex ( seq.at ( i ) ) );
        }
    }
}

int Blosum::calculateStep1Matrix ( matrix< int >& matrix1, const vector< vector<int> > & counts ) const {
    int sum = 0;
    for ( int i = 0; i < alphabetLen_; i++ ) {
        for ( int j = i; j < alphabetLen_; j++ ) {
            int result = 0;
            if ( i == j ) { // przekątna macierzy
                for ( vector< vector<int> >::const_iterator c = counts.begin(); c != counts.end(); ++c ) {
                    int cnt = c->at ( i );
                    if ( cnt > 0 ) {
                        matrix1 ( i, j ) += cnt * ( cnt - 1 );
                    }
                }
                sum += matrix1 ( i, j );
            } else {
                for ( vector< vector<int> >::const_iterator c = counts.begin(); c != counts.end(); ++c ) {
                    matrix1 ( i, j ) += c->at ( i ) * c->at ( j );
                }
                sum += matrix1 ( i, j ) * 2;
            }
        }
    }
    return sum;
}

void Blosum::calculateProbabilitiesMatrix (
    matrix< double > &      probabilities,
    vector<double> &        marginals,
    const matrix< int > &   step1_matrix,
    int                     sum
) const {
    for ( int i = 0; i < alphabetLen_; i++ ) {
        for ( int j = i; j < alphabetLen_; j++ ) {
            probabilities ( i, j ) = ( double ) step1_matrix ( i, j ) / sum;
        }
    }
    for ( int i = 0; i < alphabetLen_; i++ ) {
        for ( int j = 0; j < alphabetLen_; j++ ) {
            if ( i == j ) {
                marginals.at ( i ) += probabilities ( i, i );
            } else {
                marginals.at ( i ) += max<double> ( probabilities ( i, j ), probabilities ( j, i ) );
            }
        }
    }
}

void Blosum::calculateBlosum ( const matrix< double >& probabilities_matrix, const vector< double >& marginal_probabilities ) {
    blosum_ = matrix<int> ( alphabetLen_, alphabetLen_ );
	//zerowanie	macierzy blosum
	for (unsigned i = 0; i < blosum_.size1 (); ++ i)
        for (unsigned j = 0; j < blosum_.size2 (); ++ j)
           blosum_(i, j) = 0;

    for ( int i = 0; i < alphabetLen_; i++ ) {
        for ( int j = i; j < alphabetLen_; j++ ) {
            double fraction = ( marginal_probabilities.at ( i ) * marginal_probabilities.at ( j ) );
            if ( fraction == 0 ) {
                blosum_ ( i, j ) = 0;
            } else {
                blosum_ ( i, j ) = round ( 2 * log2 ( probabilities_matrix ( i, j ) / fraction ) );
            }
        }
    }
}


//----------------------------------------------------------------------------------------------------
//                                            TESTY
//----------------------------------------------------------------------------------------------------

void Blosum::testShowCounts ( const vector< vector<int> > & counts ) const {
    cout << "Liczba występujących symboli: " << endl;
    for ( int i = 0; i < sequenceLenght_; i++ ) {
        const vector<int> & c = counts.at ( i );
        cout << "index " << i << endl;
        for ( map<char, int>::const_iterator it = alphabetMap_.begin(); it != alphabetMap_.end(); ++it ) {
            cout << "\t" << it->first << ": " << c.at ( it->second ) << endl;
        }
    }
}

void Blosum::testShowStep1 ( const matrix< int >& substitutions, int sum ) const {
    cout << "Macierz występujących symboli: " << endl;
    cout << "\t";
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << indexToChar ( i ) << "\t";
    }
    cout << endl;
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << indexToChar ( i ) << "\t";
        for ( int j = 0; j < alphabetLen_; j++ ) {
            cout << substitutions ( i, j ) << "\t";
        }
        cout << endl;
    }
    cout << "suma: " << sum << endl;
}

void Blosum::testShowProbabilities ( const matrix< double >& m, vector<double> & marginals ) const {
    cout << "Macierz prawdopodobieństw i prawdopodobieństwa brzegowe: " << endl;
    cout << "\t";
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << indexToChar ( i ) << "\t";
    }
    cout << endl;
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << indexToChar ( i ) << "\t";
        for ( int j = 0; j < alphabetLen_; j++ ) {
            cout << m ( i, j ) << "\t";
        }
        cout << endl;
    }
    cout << endl << "\t";
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << marginals.at ( i ) << "\t";
    }
    cout << endl;
}

void Blosum::testShowBlosum() {
    cout << "Macierz BLOSUM: " << endl;
    cout << "\t";
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << indexToChar ( i ) << "\t";
    }
    cout << endl;
    for ( int i = 0; i < alphabetLen_; i++ ) {
        cout << indexToChar ( i ) << "\t";
        for ( int j = 0; j < alphabetLen_; j++ ) {
            cout << blosum_ ( i, j ) << "\t";
        }
        cout << endl;
    }
}

